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Abstract 

O 

C<^ ■ Electrical neural signalling typically takes place at the time-scale of milliseconds, and is typically modeled 

using the cable equation. This is a good approximation when ionic concentrations are expected to vary 
little during the time course of a simulation. During periods of intense neural signalling, however, the local 
extracellular K+-concentration may increase by several millimolars. Clearance of excess K+ likely depends 
partly on diffusion in the extracellular space, partly on local uptake by- and intracellular transport within 
O ' astrocytes. The processes that maintain the extracellular environment typically takes place at the time 

O . scale of seconds, and cannot be modeled accurately without accounting for the spatiotemporal variations 

in ion concentrations. This work presented here consists of two main parts: First, we developed a 
general electrodiffusive formalism for modeling ion concentration dynamics in a one-dimensional geometry, 
including both an intra- and extracellular domain. The formalism was based on the Nernst-Planck 
Cn . equations. It ensures (i) that the membrane potential and ion concentrations are in consistency, (ii) global 

particle/charge conservation, and (iii) accounts for diffusion and concentration dependent variations in 
resistivities. Second, we applied the formalism to a model of astrocytes exchanging ions with the ECS. 
r-"^ ■ Through simulations, we identified the key astrocytic mechanisms involved in K+ removal from high 

f-^ . concentration regions. We found that a local increase in extracellular K+ evoked a local depolarization 

.— iI ■ of the astrocyte membrane, which at the same time (i) increased the local astrocytic uptake of K+ 

^^ . (by locally inactivating the outward Kir-current), (ii) suppressed extracellular transport of K+, (iii) 

CO ' increased transport of K+ within astrocytes, and (iv) facilitated astrocytic relase of K+ in regions where 

the extracellular concentration was low. In summary, these mechanisms seem optimal for shielding the 
extracellular space from excess K+. 

X 

^ ; Introduction 

Electrical neural signalling is typically modeled using the cable equation, where dendrites and axons are 
represented as one-dimensional, possibly branching, electrical cables, and the transmembrane potential 
is the key dynamical variable [HE]. With the possible exception of the signalling molecule Ca^+ (see 
e.g., jH^), ion concentrations are typically assumed to be constant. The effect of ionic diffusion (due 
to concentration gradients) on the net electrical currents is neglected in standard cable theory. Also, 
resistivities (which in reality depend on ion concentrations) are assumed to be constant. These are 
often good approximations, as concentrations of the main charge carriers (K+, Na+ and CI") in the 
extracellular- (ECS) or intracellular space (ICS) typically vary little at the short time-scale relevant for 
electrical neural activity (< 100ms). 

Ion concentrations in neural tissue are actively regulated by ion pumps and membrane co- transporters. 
Regulation of extracellular ion concentrations is one of the key cellular functions of astrocytes [5^ . These 
regulatory processes typically take place at a longer time-scale (> Is) where spatiotemporal variations 
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of ionic concentrations are likely to occur. Under such circumstances, diffusion is likely to become a 
viable transport mechanism, and the standard cable- model fails to give accurate predictions . Previous 
astrocyte models have therefore been based on the Nernst-Planck equations, which consider ionic move- 
ments due to diffusion as well as electrical migration [7l[8]. It should be noted that electrodiffusion also 
may be relevant at short time-scales, e.g., in small intracellular volumes, such as presynaptic terminals 
or dendritic spines, where intracellular ion concentrations may increase by orders of magnitude within a 
millisecond [6l[9]. 

Electro diffusive models simulate the dynamics of ionic concentrations (c/c) of several ionic species (/c), 
as well as the transmembrane potential (vm)- A general concern when both vm and c^ are modeled, 
is whether the relationship between the two is consistent. Physically, vm is determined by the total 
electrical charge on the inside (or outside) of the membrane, which in turn is uniquely determined by 
the concentrations {ck) of all ionic species that are present there. In order to make specific problems 
analytically tractable, or speed up simulations, many models allow the strict dependence between Ck and 
Vm to be violated, e.g., by deriving vm from a reduced equation that does not include all ionic transports 
in the system [7^8,10 . In many situations, such approximations may be warranted. However, if applied 
to general problems, and in particular in long-term simulations, models that do not ensure an internally 
consistent Ck — vm relationship may give erroneous predictions. For example, in some heart cell models, 
ion concentrations have been reported to drift to unrealistic values in long-term simulations, while vm 
maintain realistic values pJHT3] . 

Qian and Sejnowski developed a model for electrodiffusion in excitable cells, which in a consistent 
way derives vm from the intracellular ionic concentrations [6 . Like the standard cable model, the 
electro diffusive model assumes that transport phenomena are essentially one-dimensional. Unlike the 
standard cable model, the electrodiffusive model includes the concentration dynamics of all involved ions, 
the diffusive currents arising from intracellular concentration gradients, and the concentration dependent 
variation of the intracellular resistivities. However, an important limitation with this model is that it 
only includes intracellular transports, whereas the ECS is assumed to be isopotential and with constant 
ion concentrations. In reality, the ECS comprises about 20% of the total neural tissue volume, while the 
remaining 80% is the ICS of various cells (mainly astrocytes and neurons). When a large number of cells 
participate in simultaneous ion exchange with the ECS, the impact on the ion concentrations in the ICS 
and ECS may be of the same order of magnitude. For example, it is known that the local extracellular 
K+-concentration can increase by several millimolars (i.e., more than double) during periods of intense 
neural activity [7l[l4l[T5]. Furthermore, clearance of excess K+ from high concentration regions likely 
depends partly on diffusion in the ECS, partly on local uptake via astrocytic K+-uptake mechanisms, 
and partly by intracellular transport within astrocytes [7|[8| [T5i rr6] . An understanding of these processes, 
requires a general electrodiffusive framework that explicitly includes both the ECS and the ICS. 

In situations when macroscopic transport processes are effectively one-dimensional, the complex com- 
position of the tissue (Fig.[T]A) can be simplified to the two-domain model shown in Fig.[T]S [7l[8]. There, 
the ICS of all cells participating in the transport process have been represented as an equivalent cable 
(/-domain), which is coated by ECS (£^-domain). The I-E system may be pictured phenomenologically 
as an average single cell coated with the average proportion of available ECS per cell. Such a geometri- 
cal simplification was previously motivated for one-dimensional transport phenomena through the glial 
syncytium [ZlIH]. 

In this work, we first derive a simple, general mathematical framework for modeling the dynamics of 
the membrane potential (vm)^ the intra- (cki) and extracellular (ckE) ion concentrations for a set {k) 
of ionic species, and identify the conditions where the formalism reduces to the standard cable model. 
Next, we apply the electrodiffusive formalism to model ionic exchange between astrocytes and the ECS, 
and investigate the relative role of astrocytes in K+ removal from high concentration regions. Finally, 
we provide a discussion of our results, and of the underlying assumptions in our new electrodiffusive 
framework. 




Figure 1. A two domain- model for ion concentration dynamics in the intra- and 
extracellular space, when macroscopic transport is essentially one-dimensional. (A) A piece 
of neural tissue with cross section area Aref and an arbitrary extension / in the x-direction. The tissue 
contains cehs (dark grey) that participate in the transport process, and cells that do not (light grey). , 
where aj is the fraction of Aref that is the ICS of the participatory cell type, qe is the fraction that is 
ECS, (B) The participatory cells represented as an equivalent cylindrical cable (/), coated by ECS (E). 
The geometry is specified by three parameters, where aj and aE are, respectively, the fractions of Aref 
occupied by the ICS of participatory cells and the ECS, and OM{'m~^) is the amount of membrane area 
per tissue volume, or, equivalent ly, the circumference of the equivalent cable. Due to the presence of 
non-participatory cells, we generally have that aj -\- aE < 1. 

Results 

Electrodiffusive formalism 

A formalism was derived for computing the ion-concentration dynamics in a geometry as that depicted in 
Fig.[T]B. The formalism is summarized in Fig. [21 and the derivation is included in the following subsections. 



Particle conservation 

In Fig. [T]B, particles in / or £^ may move along the x-axis or across the membrane. In a segment Ax 
of /, centered at x, and with volume a/ Ax, the particle concentration dynamics of an ion species k is 
determined by: 



-OM^xjkM{x,t) -^aijkiix- Ax/2,t)) - aijki{x ^ Ax/2, t) = ajAx ^^ ' , (1) 

where the transmembrane- {jku)-, the intracehular- {jki) and the extracehular {Jue) flux densities of 
particle species /c, have units mol/(m^s). The flrst term on the left represents the ionic flux that enter (+) 
this segment through the piece of the membrane with area OmAx. The second and third terms represent 
the ionic fluxes that enter (+) /leave (—) the section through the left/right boundaries, with cross section 
areas a/. If the net flux into the segment is nonzero, the ion concentration will build up over time, 
according to the right hand side of Eq. [TJ 

We divide Eq. [T]by ajAx, and take the limit Ax -^ 0, to obtain the continuity equation on differential 
form: 

djki{x,t) Om . / ,x , dcki{x,t) 

^d^ + —JuM{x,t) + ^^^ = (2) 

^^^-^..M(.,t) + ^^^ (3) 

ox aE ot 

We have also written up the continuity equation for the extracellular domain. By convention, jku has 
been defined as positive in the direction from / to E. 

The axial flux densities are described by the generalized Nernst-Planck equation: 

. Dkdckn{x,t) DkZk . dvn{x,t) 



where Zk is the valence of ion species /c, and the index n represents I or E. The flrst term on the right 
in Eq. |4]is the diffusive flux density (j^^), driven by the concentration gradients, and the last term is the 
fleld flux density (j^^), i.e., the flux density due to ionic migration in the electrical fleld. The effective 
diffusion constant Dl = Dk/X^ ^^ composed of the diffusion constant D^ in dilute solutions and the 
tortuosity factor An, which summarizes the hindrance imposed by the cellular structures 0117]. We 
use i/j = RT/F{mV), where R = 8. 3144621 J/(mo/iir) is the gas constant, T the absolute temperature, 
and F = 96, 4853365(7/77io/ is Faraday's constant. As in [8], we have then implicitly assumed that the 
Einstein relation holds between the effective diffusion constant and effective electrical mobility. 

The formalism is general to the form of jkM •, which may include contribution from multiple membrane 
mechanisms, such as ion pumps, co- transporters and ion channels. It is sufficient to require that jkM is 
known at any point in time given the voltage across the membrane, the ionic concentrations on either 
side of the membrane, and possibly some additional local information (mi,m2,etc.) reflecting the local 
state of the membrane: 

jkM{x, t) = f{cki{x, t),CkE{x, t),VM{x, t), fai (x, t), m2(x, t), ...). (5) 

As boundary conditions, we shall apply the sealed-end condition, i.e., we assume that no fluxes enter 
or leave through the ends (x = and x = I) oi I or E\ 

jkn{0,t)=jkn{l.t)=0. (6) 

Equations [2] - [31 together with with Eqs.|ll[5]and [6l specify the system we want to solve. Before we 
derive the electrodiffusive formalism for this problem, we recall how the standard cable equation can be 
derived from the principles of particle conservation. 



Charge conservation 

The particle conservation laws (Eqs. [2]- [3]) can be transformed to charge conservation laws by the use of 
the general relations (see e.g., [18 J: 

Pn{x, t) = F^ ZkCkn{x, t) + psn{x) (7) 

k 

iM{x,t) = F^ZkjkMJx.t) (8) 

k 

in{x,t) = F^ Zkjkn {X,t). (9) 

k 

Here, pn{C/m.^) is the charge density, Im (A/m^) is the transmembrane current density, and in (A/m^) 
is the axial current density. For practical purposes, we have included a density of static charges {psn) in 
Eq. [71 representing contributions from ions/charged molecules that are not considered in the conservation 
equations. If the set Ckn include all present species of ions, then psn = 0. To keep notation compact, we 
from here on omit the functional arguments (x,t). 

If we multiply the particle conservation laws (Eqs. [2]- [3]) by Fzk^ take the sum over all ion species, /c, 
and use Eqs. [71 -[H we obtain the equivalent laws for charge conservation: 

ai^ + OMiM + ai^ = (10) 

dE-^ Om^M + ^E-^ = 0. (11) 

Note that the last term only depends on the mobile ions, as dpsn/dt = 0. 

Standard cable equation 

The standard cable equation may be derived by combining the charge conservation laws (Eqs. [TQ1-[TT]) 
with three simplifying assumptions: (i) E is assumed to be isopotential and with zero resistivity, (ii) the 
membrane is a parallel-plate capacitor, and (iii) ion concentrations are effectively constant, i.e., diffusive 
currents are negligible and resistivities (see Eq. [T5l below) are constant. 

Assumption (i) implies that we only need to consider charge conservation in / explicitly. To obtain 
the cable equation in the standard form, we must express pi and ij in Eq.[TO]in terms of '^m and Ovm/Ox. 

Assumption (ii) allows us substitute vm for pj. A capacitor with capacitance SC separates a charge 
SQ from the opposite charge —SQ^ and generates a voltage difference v = SQ/SC. The charge inside a 
piece {Sx) of membrane with area Om^x is SQj = pjajSx. The capacitance of this piece of membrane is 
SC = CmOmSx^ where Cm denotes the membrane capacitance per membrane area. We therefore obtain: 

_ SQi _ piajSx _ _a^_p^ 
SC CmOmSx Dm Cm 

According to assumption (iii), diffusive currents are negligible, and Eq. [H reduces to: 

./ DkZk dvi 

J/e/=J.. = -^^c.. — . (13) 

If we insert Eq. [13] into Eq. [9l we see that the axial current density obeys Ohm's current law: 
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FDkzl dvi 1 dvi 



^ Xjip dx rj dx ' 



where we have identified the resistivity, rni^m): 



in the ICS {n = I). Finahy, we insert Eqs. [T4l and [T2l into Eq. [TO]to obtain the cable equation: 

-^^M ^Cm^z— =0. (16) 



Om ^i dx'^ dt 

Note that r^ generahy depends on Ckn- However, we have here assumed that r^ is constant (cf. assumption 
(iii)). Furthermore, we have used the identity: dvj/dx = Ovm/Ox, which fohows from the definition 

vm = vi -ve, (17) 

together with the assumption (i) that E is isopotential. Eq. [16] is the most commonly used form of the 
cable equation, although there are versions that also explicitly considers the extracellular domain [2] . 

Two-domain electrodiffusive model 

The cable equation only considers the net electrical transports, and "hides" the underlying transports of 
different ionic species. We now develop the electrodiffusive formalism for computing the ion-concentration 
dynamics. Like in standard cable theory, we limit the study to the one-dimensional geometry in Fig. [T]B. 
Unlike standard cable theory, we explicitly consider both domains / and E^ and we do not neglect diffusive 
currents nor concentration dependent variations of the resistivities. 

The conservation equations (Eqs. [2] -[3]), with the Nernst-Planck equation (Eq. H]) for jkn specify the 
system we want to solve. As in standard cable theory, the formalism is general to the form of jkM (Eq.[5]). 
With TV ion species, Eqs. [2] - [3] represent a system of 2N + 3 variables which are functions of x and t. 
These are the 2N concentration variables {ckn for /c = 1,2, ...N and n = E^I)^ and the three additional 
variables {vu^dvi /dx and dvE/dx) occurring in the expressions for the fiux densities. 

To reduce the number of independent variables to the 2N state variables {ckn) we need three conditions 
relating vm^ dvi/dx and Ove/Ox to Ckn- The first two conditions we recognize from standard cable theory, 
namely (CI) that vm is determined by the charge density (Eq. [12]), and (C2) the general definition of 
(Eq. [TT]) of vM' In the two-domain case, we use the additional condition (c3): 

aipi = -aEpE' (18) 

We shall refer to C3 as the charge symmetry condition. Its origin is explained below. 
According to condition CI, vm is given by Eq. [T2l 

where we have inserted Eq.[7]for p/, so that vm is expressed in terms of ionic concentrations. Equivalent ly, 
we may also express vm in terms of the ion concentrations in the ECS: 

VM = - ^ f. PE = - ^ ^ {F'^ZkCkE^PsE), (20) 

where the negative sign follows from the convention that vm is positive when / is positively charged. By 
demanding consistency between Eq. [19] and Eq. [20] we can derive Eq. [TH] which is the charge symmetry 
condition (C3). It implies that the charge on the inside of a piece of membrane is equal in magnitude and 



opposite in sign to the charge on the outside. CI and C3 are both imphcit when the membrane is assumed 
to be a parallel plate capacitor. C3 is also related to the issue of electroneutrality (see Discussion). 

The next step is to express the voltage gradients (dvn/dx) in terms of ionic concentrations. The 
constraints C2 (Eq. [T7|) and C3 (Eg. [T8|) allow us to derive two independent equations that relate dvE/dx 
and dvj/dx. The first equation is obtained by differentiating Eq. [171 

dvMJx) ^ dvi{x) _ dvEJx) 
dx dx dx 

We recall that vm is already a known function of ion concentrations (Eq. [19] or Eq. [20]) . 

A second equation relating dvj/dx to Sve/Ox may be derived by combining Eq. [18] with the charge 
conservation laws. If we sum Eqs.[TQ]and [TT| we immediately see that the terms involving im cancel out. 
Due to Eq. [THl also the last terms on the left cancel, so that we are left with: 

Due to sealed end-condition (Eq. [6]), in(0) =0, so that Eq . [22] takes the simple form: 

aiii = -aEiE- (23) 

If the charge symmetry condition (C3) is satisfied at a given time t = (and we must specify the initial 
concentrations so that this is true), Eq. [23]is the condition that it is satisfied at all times t. 

We now decompose the current density in a diffusive term and a field term: in = in-\-il, and express 
z{ in terms of Ohm's law (cf. Eq.[14]). If we insert this into Eq.[23l we obtain the second equation relating 
Ove/Ox and dvi/dx: 

■d , 1 dvi\ _ _ (^^ , I dvE 



Finally, Eq. [21] and Eq. [24] can be solved for the voltage gradients. After some simple algebra we 
obtain: 

ox \ ox aE J \ riaE J 

dvE f dvM , . -d , riaE ,d\ ( , , riaE\~ 



ox \ ox cii J \ rEdi J 

Here, r^ is given by Eq. [151 ^n t>y Eq. [H and vm by Eq. [19] or Eq. [20] All voltage terms are thereby 
expressed in terms of ionic concentrations. With this, the conservation equations (Eqs. [2] - [3]) are fully 
specified, and can be solved numerically with appropriate boundary conditions. The final set of equations 
is summarized in Fig. [2] 

Electrodiffusive formalism vs cable equation 

From Eq. [TO]for charge conservation in /, we may derive a differential equation for the time development 
of vM' We use Eq. [T9]to substitute vm for pi. Furthermore, we use the decomposition ij = ij -\-i^j^ with 
Eq.[llfor i{, and Eq.[25]for ^^^. We then obtain: 



dx 

aj d 



ri ri dx J \ riaE J 



Ij Ie 



»m + Cm^=0. (27) 

ot 



Om dx 
This is the equivalent to the standard cable equation (Eq.[T6|). for the electrodiffusive two-domain system 



Specify initial conditions: ^mo and Ckno for n = I,E and /c = 1, 2, 3, . . . , A^ 
Define static cliarge densities: 

PsI = C'm^mO — F{ckIO + CiVa/O " Cc^/q) 

a/ 

PsE = CmVmO " F{ckEO + CAra£;0 " CciEo) 

dE 

Solve the 2N conservation equations: 

djki{x,t) Om . / ,x , dcki(x,t) 
-^^ + ^Jkuix, t) + -^^ = 



with: 



. Dkdckn{x,t) DkZk , dvn{x,t) 



^2 



^ = E^^^^^ 



j;i(a:) = -F^ 



fc 



ZkPk dCkn{x,t) 

\l dx 



vm{x)= i^ {Fy2^kCki + Psi) or - n n i^yZ^kCkE + Pse) 



Figure 2. The two-domain electrodiffusive formalism. The set of equations summarizes (and 
fully specify) the electrodiffusive formalism. It is applicable to general problems. The transmembrane 
currents (j/cm) need to be specified for any membrane mechanism included in a model. 

A few notes: Firstly, a corresponding dynamical equation for vm could have been derived from the 
extracellular conservation law (Eq.fTT]). Due to the charge symmetry condition, the two equations would 
be equivalent. Secondly, unlike the standard cable equation, Eq. [271 does not provide a complete system 
description, as Eqs.[2]-[3]must be solved to determine i^ and r^. Thirdly, when the ionic concentrations 
are known, Eq. [271is not necessary for computing vm^ as vm can be computed algebraically from Eq.fTOl 
Eq. [23 is mainly useful for comparison with the standard cable equation. 

We can immediately see that if we make the common assumptions (i) that the extracellular resistivity 
{ve) is zero, (ii) that the diffusive currents (i^) are zero, and (iii) that the intracellular resistivity (r/) 
is constant, then Eq. [27] reduces to the standard cable equation (Eq. [T6|). We should note that there are 
two-domain versions of the cable equation where the first assumption is not made [2]. The two other 
assumptions are warranted only in cases when the spatiotemporal variations in ionic concentrations is 
such that r/ varies little, and i^ <C i^ during the time course of a simulation. 



Astrocyte Model 

The electro diffusion- formalism was applied in a one-dimensional model for astrocytes exchanging ions 
with the ECS. The purpose was to investigate the relative role of astrocytes in K+ removal from high 
concentration regions. 

The model was developed for macroscopic transport processes, involving all astrocytes in a piece 
of tissue. The geometry in Fig. [T]B was therefore applicable, with / representing a phenomenological 
"average" astrocyte (the cable, /), surrounded by a sheet of ECS (the coating, E). The geometrical 
parameters a/, aE and Om have been estimated for astrocytes in neural tissue (see Table 1). For the 
extension in the x-direction, we used / = 300/im. 

Astrocytic membrane mechanisms were adopted from a previous point-model of an astrocyte [16] . The 
included mechanisms were the standard, passive Na+ and CI" channels, the inward rectifying K+-channel 
(Kir), and the Na+/K+-pump, as sketched in Fig.[3l The membrane mechanisms are described in further 
detail in the Me ^/loc/^- sect ion. 
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Figure 3. Astrocyte model. A representative astrocyte (/) exchanging ions with the ECS {E) 
through membrane mechanisms as indicated. The system input to the ECS was applied in the input 
zone. The output was applied over the whole axis. The decay zone was defined as the part of the axis 
where no input was applied. 



We assumed that only the three main charge carriers (K+, Na+ and CI") contributed to electrodiffusive 
transport. For the diffusion constants (I^/c), we used values valid for electrodiffusion in diluted media [19], 
modified with the tortuosities (A^) estimated in [8]. The same values have also been used in earlier, related 
studies [6l[20]. All relevant model parameters are listed in Table 1. 

As external input to the system (mimicking the effect of enhanced, local neuronal activity), we used 
a constant K+/Na+-exchange (adding K+, subtracting Na+) in a selected region (0 < x < I/IO) of the 



10 



ECS between t = 100s and t = 400s: 

JK = -JNa = Jin iov < X < I/IO & 100 < t < 400 



(28) 



Conversely, the model output was a concentration- dependent Na+/K+-excliange (subtracting K+, adding 
Na+), occurring over the entire astrocyte axis, causing the extracellular K+-concentration to decay to- 
wards the resting concentration. The output could represent the uptake of K+ /release of Na+ by other 
neurons via K+/Na+-exchangers and other uptake mechanisms: 

jout ^ _jout ^ -kdec{cK - CKo) for ah X & t .^g) 

The decay factor {kdec) was set to a realistic value for maximal neuronal K+/Na+-exchange under phys- 
iological conditions (see Table 1). The input flux density {jin) was specified to a value that gave a 
steady-state K+-concentration of about lOmM in the input region {0 < x < I/IO) during constant input 
condition (see Fig. H]). 

Ion concentration dynamics in the Astrocyte/ECS system 

Fig. mA-E" shows the dynamics of the atrocyte/ECS system in the middle of the input zone {x = 1/20). 
During the input {100s < t < AOOs) there was a net influx of K+ and a net efllux of Na+ to/from the 
input zone (Fig.|4]A). The constant input caused an increase in ck and a decrease in CNa both in / and E 
(Fig. |4]B-(7). Although CI" was not added to the system, ecu increased on behalf of ccie- The changes 
in ionic concentrations coincided with a depolarization of the membrane in the input zone (Fig. |4]E), 
reflecting concentration dependent changes in the reversal potentials of the involved ionic species. 

The total output rate (integrated over all x) increased with time, due to the general increase in 
ckE' When cke became sufficiently high, the total output rate coincided with the total input rate, and 
the system reached steady-state (SS). It took 49 s from the constant input had been turned on until the 
slowest variable {\Accie\) reached 99% of its saturation value. Most variables approached SS signiflcantly 
faster (e.g., 12 s for \Acke\ and 19 s for \Avm\)' When the input was turned off, the system gradually 
returned to the original resting state. 

We here focus on the SS-situation, i.e., on the activity of astrocytes during periods of on-going intense 
neural activity. Fig. |4]F- J shows the spatial profiles of the input/output, the ionic concentrations and 
the membrane potential at a time tss = 400s, when the system was in SS. For all variables, deviations 
from the resting values were greatest in the input zone. There, cke was about 10 mM during SS, i.e., 7 
mM above the resting concentration (3 mM). 

The resting potential vmo = — 83.6m V corresponded to concentrations Cei — O.lSmM and Cqe — 
0.36mM of unit charges in the ICS and ECS (cf. Eq. [39]). At SS, vm had increased from the resting po- 
tential to about -60 mV, consistent with small absolute changes (Acg/ ^ 0.05mM and AceE — O.lOmM) 
in the concentration of unit charges. As seen in Fig. |4]B-(7, these changes were very small compared to 
the changes Ac^n ^^ any of the ionic concentrations. Variations in ion concentrations were thus always 
so that anions and cations remained closely balanced in numbers, giving rise to a relatively small net 
charge. 

The concentration dependent changes in the resistivities were quite significant. In the input zone, r/ 
decreased by about 10%, and te increased by about 20% compared to their basal values (Fig. |4]E). 

Ion transport pattern in steady state 

During SS, the system output was distributed over the x-axis, with equal areas under the input and 
output curves (Fig. |4]F). In the input zone, the output rate was about 1/3 of the input rate. This means 



Table 1. Model parameters 
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Parameter 

/ (length of astrocyte) 

ai (astrocyte volume/total tissue volume) 

aE (ECS volume/total tissue volume) 

Om (glia membrane area/total tissue volume) 

Dk (K+ diffusion constant) 

DNa (Na+ diffusion constant) 

Dci (CI" diffusion constant) 

A/ (intracellular tortuosity) 

A^; (intracellular tortuosity) 

Cm (specific membrane capacitance) 

Qkq (baseline K+-conductance) 

QNaQ (baseline Na+-conductance) 

gciQ (baseline Cl"-conductance) 

Pmax (maximum Na+/K+ pump-rate) 

Kke (c/c£; -threshold for Na+/K+ pump) 

KNai (cATaJ-threshold for Na+/K+ pump) 

cJ^^Q (initial ECS K+-concentration) 

^NaEO 



(initial ECS Na+-concentration) 
(initial ECS K+-concentration) 



'-Naio (iiiitial ICS Na+-concentration) 
^ciEO (initial ICS Cl"-concentration) 
c^ijQ (initial ICS Cl"-concentration) 
v'^Q (initial membrane potential) 
^dec (decay factor for cke) 
jin (constant input in input zone) 



Value 

300/i771 

0.4 
0.2 

4.8 X lO^m-i 

1.96 X 10-^m^/s 

1.33 X 10-^mVs 

2.03 X 10-^mVs 

3.2 

1.6 

ljj.F/cm^ 

O.bS/m'^ 

1.12 X 10-^mol/{Sm^) 

l.bmol/m? 

lOmol/m^ 

3.0 + 0.082mM 

100.0 - 0.041mM 

145.0 - 0.378mM 

15.0 + 0.189mM 

134.0 - 0.29mM 

5.0 + 0.145mM 

-85 + lAmV 

2.9 X lO-^m/5 

7 X 10~^mol/{m'^s) 



Reference 

E] 

[6!fT9!l20l 
[6,19,20^ 
[611191120] 
i 

la 
m 

M 
M 
M 

la 
la 
ca 
la 
ca 
ca 
ca 
ca 
ca 



Initial concentrations are given as CknO = Value from flE\ + Correction^ where the sum gives rise to a 
system at rest. 

^ The maximum average Va+/i^+-pump rate for a single neuron was estimated to 
A = 2 X 10~ ^ mol / {m^ s) [21j. We obtained kdec by solving kdec{c^^E ~ ^keo) = ^^ assuming that 



-KE 



= lOmM. 



that about 2/3 of the K+ that entered the system was transported in the positive x-direction, and left 
the system from the decay-zone (cf. Fig. [3]). 

We wanted to explore the role of the astrocyte in removing K+ from the input zone, relative to the 
role of axial transports in the ECS. To do this, we analyzed the spatial profiles of all ionic flux densities 
during SS (Fig. [5j). We distinguished between field flux densities (j^^) and diffusive flux densities (j^^). 
The main transport routes during SS are summarized Fig. [5]G^. The net CI" transport {Jcie ~^ Jcie) ^^^ 
very small, and was not included in the summary. 

For Na+ and K+, the main transport routes were as follows: Na+ entered the system in the decay 
zone of the ECS, was transported into the input zone, and left the system from the input zone. The 
main axial transport occurred in the ECS. In contrast, K+ entered the system in the input zone, where 
a major fraction of it crossed the membrane. Transport of K+ out from the input zone predominantly 
took place inside the astrocyte. Outside the input zone (i.e., in the decay zone), the astrocyte released 
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Figure 4. Dynamics and steady state profiles for the astrocyte/ECS-system. (A-E) 

Dynamics of selected variables in a point (x = 1/20 = 0.15/im) in the middle of the input zone. The 
input was given from t = 100s to t = 4005. (F-J) Spatial profiles of selected variables at a time 
{t = 400s) when the system was in steady state. (B-C, G-H) Ionic concentrations are represented in 
terms of deviations from resting concentrations: Ackn = Ckn — CknO for n = I^E. (E,J) Resistivities are 
plotted relative to the resting values rEo = lAb^m and r/o = 12.0rtm. 



K+ to the ECS, from where it eventually left the system. The sum of the Na+ and K+ transports gave 
rise to a net electrical current which cycled in the system. 

Two basic mechanisms explain the qualitative difference between Na+ and K+ transports. Both 
are related to the membrane being most depolarized in the input zone (Fig. |4]F). The first mechanism 
concerns the transmembrane fluxes. The Na+/K+-exchanger mediated an inward flux of K^ and an 
outward flux of Na'^ . The exchanger was counterbalanced by passive fluxes in the opposite direction 
{Na~^ in and K^ out), proportional to vm — e^. In the case of A^a+, the passive flux and the exchanger 
rate were closely balanced across the length of the astrocyte, so that JNaM was small everywhere. In the 
case of i^+, vm — ^k was very small in the input zone (e^ ^ —62mV (Eq. [34]) and vm ^ — 60mF), but 
quite big in the decay zone. In the input zone, Jkm was therefore dominated by K+-uptake through the 
Na+/K+-exchanger, while i^+-release through the outward Kir-channel dominated in the decay zone. A 
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similar role of the Kir-channel in spatial K+-buffering has been suggested previously [22^. 

A second (and to our knowledge, novel) mechanism for explaining the differences between the Na+ and 
K+ transports concerns the axial fluxes. As the astrocyte was most depolarized in the input zone, / had 
the highest positive charge density and E had the highest negative charge density there. Therefore, the 
electrical force on Na+ and K+ (being cations) was in the negative x-direction in the ECS {dVE/dx > 0), 
and in the positive x-direction in the ICS {dVj/dx < 0). Transport of K+ out from the input zone 
therefore had the best conditions in /, where diffusive and electrical forces were driving the ions in the 
same direction. For the same reason, Na^ transport into the input zone dominated in E. 

In summary, the local depolarization of the astrocyte induced changes in both the transmembrane 
and axial flux densities, which both improved the efficiency of the astrocyte in removing K+ from the 
input zone. 

Consistency of formalism 

In the simulations above, vm was defined in terms of the charge density in /, and computed algebraically 
by solving Eq.[T9]at each time step. Identical results (down to a very small numerical error) were obtained 
when Vm was defined by the charge density in E (Eq. [2Q|) . and when vm was computed differentially by 
using Eq. [23 (results not shown) . As all transports are included in Eq. [271 the algebraic and differential 
methods yielded consistent results. 

When the input was turned on and off, a small numerical error was introduced in the conservation 
of ionic concentrations, inducing a small error in the total charge in the system. The relative deviation 
from global charge neutrality {Qtot = 0), defined as CQtot = Qtot/{\QE\ + |Q/|), where Qj and Qe refer to 
the total charge in / and E^ was about 10~^^. This gave rise to a relative deviation from perfect charge 
symmetry (cf. Eq. [TH]), defined as Cp = {ajpi — aEpE)/(\ciipi\ + \(IePe\)^ which was also on the order of 
10~^^ (for all x). Accordingly, vm computed from the charge density in E deviated by a relative factor 
~ 10~^^ from Vm computed from the charge density in /. This corresponded to an absolute difference 
of ~ 10~^7TiV. Errors tended to be somewhat larger when the differential method was used. Then vm 
deviated locally by up to ~ l{)~^mV from vm derived from the charge density in / or E. 

Errors will generally depend strongly on the algorithm used for solving the differential equations, the 
time step, and the number of compartments in the simulated system. The errors could likely be reduced 
by using a smoother input signal than the step function in Eq. [28l We did not engage in further analysis 
of the origin of the errors, but were content with their smallness. 

Discussion 

We presented a one-dimensional, electrodiffusive framework for modeling the dynamics of the membrane 
potential {vm) and the ion concentrations {ckn) of all included ion species [k) in an intra- and extracellular 
domain (Fig. [2]). The framework could have a broad range of applications within the field of computational 
neuroscience. In the current work, it was applied to simulate the role of astrocytes in K+-removal from 
high concentration regions. 

Spatial K+-buffering by astrocytes 

The astrocyte/ECS-model provided a mechanistic understanding of how astrocytes may remove K+ from 
high-concentration regions. In summary, the model astrocyte responded to a local extracellular increase 
in ck by a local depolarization of the membrane. At the same time, this depolarization (i) increased 
astrocytic K+ uptake in the input zone, (ii) increased astrocytic K+-release outside the input zone, (iii) 
decreased axial K+ transport in the ECS, and (iv) increased axial K+ transport inside the astrocyte. 
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Figure 5. Transports in the astrocyte/ECS system during steady state (t =400s). (A-B) 

Flux densities in E and / due to electrical field. (C-D) Flux densities in E and / due to diffusive forces. 
(A-D) Flux densities jkn were scaled by the relative area fraction a^ so that fluxes in / and E could be 
compared directly (I.e., aEJkE = CLjjki, gives the same total flux of k in I and E). (E) Total flux 
densities into system {input — output). (F) Transmembrane flux densities. (G) Flow chart showing the 
main transports during steady state. The length of the arrows indicate flux densities, but are not 
numerically exact. 



Above, (i-ii) directly concern well documented astrocytic membrane mechanisms. Increased uptake 
in high concentration regions (i.e., the input zone in Fig. [3]) was effectively achieved by inactivation 
of the outward Kir-channel, facilitating a high net uptake through the K+/Na+-exchanger. This is 
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in agreement with previous experimental findings [5l[23]. In addition (iii-iv), we also found that the 
astrocyte induced changes in the intra- and extracellular voltage gradients that facilitated intracellular 
K+-transport. Regulation of the longitudinal transport represents a (to our knowledge) novel mechanism 
that astrocytes may utilize to shield the extracellular space from excess K+. All these effects (i-iv) taken 
together turned the astrocyte into an efficient sluice for removing K+ from the input zone. 

The input to the astrocyte/ECS model was a K+/Na+ exchange in the ECS. Within a few seconds, this 
evoked large concentration gradients in the ECS and ICS. These circumstances differ from the constant- 
background conditions assumed in most neural models, but are physiologically realistic under periods of 
intense neural activity [8^14 . As we saw in Fig. HI ionic concentrations changed by several mM before the 
system reached SS, while Cqe and Cg/ remained bounded by the extreme electrical forces that would be 
associated with high absolute values of vm- When the system reached SS, the resistivity had increased by 
up to ~ 20% from the resting value, and the diffusive current in the ECS was about 25-30% of the field 
current. Hence, the predictions made by the electrodiffusive model differed significantly from what would 
be predicted by standard cable theory, where diffusive currents and concentration dependent variations 
in resistivities are neglected. 

Astrocytes are known to possess several membrane mechanisms that were not included in the current 
model. K+-uptake by Na — K — 2C/-cotransporters and K — C/-cotransporters [16^ are two candidate 
mechanisms that could affect the simulated results. 

Macroscopic transports vs. single cell models 

The astrocyte/ECS-model was represented phenomenologically as a single astrocyte coated with the av- 
erage proportion of available ECS per astrocyte (Fig. [3]). This geometrical representation is motivated 
for macroscopic transport processes, when a large number of astrocytes perform the same function simul- 
taneously [8]. For spatial K+-buffering, this was a reasonable assumption, as the input was a change in 
the ion-concentrations in the ECS, shared by all present astrocytes. 

If we instead wanted to study a cell specific signal, such as the response of a single astrocyte to a 
transmembrane current injection, the geometrical representation in Fig. [T]B would be less appropriate. 
Firstly, the notion of the ECS as a relatively thin coating following a single cell is only motivated at 
the macroscopic "average transport" -level. Secondly, a/ and Om would refer to the single participatory 
astrocyte, and not all astrocytes intersected by Aref- Then we would expect that a/ << a^;, as a single 
active cell would have a significantly larger proportion of the ECS to its own disposal. In single-cell 
models it is common to assume that conditions in E are constant, so that only / is modeled explicitly. In 
this limit, the electrodiffusive formalism reduces to the one-domain model presented previously by Qian 
and Sejnowski 0. 

Relationship to other electrodiffusive modeling schemes 

The framework presented here is essentially an expansion of the one-domain model by Qian and Sejnowski 
[6] to a two-domain model that includes both the ECS and ICS. Like the one-domain model, the framework 
ensures (i) a consistent relationship between vm and Ck- Unlike the one-domain model, the framework 
ensures (ii) global particle/charge conservation, and (iii) that the charges on either side of a piece of 
membrane must be equal in magnitude and opposite in sign {5Qi = —SQe)- The latter constraint is 
implicit when the the membrane is assumed to be a parallel plate capacitor, an assumption made in most 
models of excitable cells (see e.g., [ll[2l[6l[18] ) . It is also related to the topic of electroneutrality. 

Electroneutrality in electrodiffusive models of biological tissue has been the topic of many discussions 
[24H26] . It is relevant for how the electrical potential {v)^ occurring in the Nernst-Planck equation, is 
derived. Generally (at sufficiently course spatial resolutions so that the charge density can be assumed 
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to be continuous), v obeys Poisson's equation: 

V ■ (eV^;) = -p, (30) 

where e is the dielectric constant, and p = ps -\- F^^{zkCk) is the total charge density. 

In biological tissue, the charge relaxation time t = Vn^n is very small in any region except in the thin 
Debye layer (~lnm) surrounding a bio- membrane. Any nonzero p in the bulk solution will decay very 
rapidly (r ~lns) to zero [19]. Several models have simulated electrodiffusion by solving the Nernst-Planck 
equations in one or more dimensions, with Poisson's equation for v (see e.g., [25 1I27H3Q] ) . The advantage 
with this procedure is that the Poisson-Nernst-Planck (PNP) equations can be implemented in a general 
way in three-dimensional space. The challenge is then to specify the appropriate boundary conditions 
for solving Eq. [30] in the vicinity of membranes. Generally, PNP-solvers apply a fine spatial resolution 
near the membrane, and simulation time steps smaller than the charge-relaxation time ^2^. For these 
reasons, they tend to be extremely computationally demanding [3T] . 

The formalism presented in this work belongs to a class of of one-dimensional models, including the 
cable equation and several electrodiffusive models [6-8,10,32 , which bypasses the computationally heavy 
PNP-scheme. The physical interpretation of these models is as follows: Any net charge in a volume 
An^x is implicitly assumed to be located in the thin Debye-layer surrounding the capacitive membrane, 
and is identical to the charge that determines vm- The remainder of the space (i.e., the bulk) will 
therefore be electroneural {ptot = 0). Note that any finite volume, enclosing a piece of membrane, will 
also be electroneutral. This follows from the charge symmetry condition (Eq. [T8|) . constraining the charge 
on either side of the membrane to be equal in magnitude and opposite in sign. The charge symmetry 
condition and the electroneutrality condition are in this way closely related. In these electroneutral 
models, charge relaxation is implicit. This is a plausible assumption at time scales relevant for most 
biophysical processes. Accordingly, simulations may be run with time-steps ranging from 1 ms to 1 s, 
depending on the time course of the included membrane mechanisms. 

To our knowledge, the formalism summarized in Fig. [2] is the first biodiffusive model where the intra- 
and extracellular voltage gradients have been derived from the charge symmetry condition. Eqs. [25land 
[261 can be interpreted as summarizing all local and global electrical forces driving the system towards 
electroneutrality. 

A natural future ambition would be to generalize the electrodiffusive formalism to 2 or 3 spatial 
dimensions, so it can address the same 3-dimensional transport problems as PNP-solvers. The challenge 
will be to formulate the system as a grid of coupled constraints (electroneutrality in the bulk and EqJT2l 
for Vm across the membrane) for which the Nernst Planck-equations can be solved with time steps much 
longer than those involved in the charge relaxation process. 

Model and Methods 

Astrocytic membrane mechanisms 

The transmembrane ion fluxes in the astrocyte model were: 

QKQJKir / \ oo /o1^ 

JKM = [vm - ex) - 2P (31) 

r 

JNaM = —^ {vm - eNa) + 3P (32) 

r 

jciM = ^^ {vm - eci) ' (33) 

r 

Here, g^o are the baseline conductances of the passive K+, Na+ and CI" currents. The currents depend 
linearly on the difference between vm and the reversal potential, 

e/e = ^log(c/e£;/c/e/), (34) 



for the respective ion types (k). The potassium current was modified by the Kir- function |8]: 

l + exp(18.4/42.4) 



fKir{cKE,^V,VM) = xl 

ckeq 



l + exp[-(118.6 + eKQ)/44.1] 
l + exp[-(118.6 + ^M)/44.1] 
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(35) 



_l + exp[(Av + 18.5)/42.5]_ 

where A'^ = Vm — ^k^ and exo(^^) is the Nernst potential for potassium at basal concentrations CkEQ 
and c/e/o. 

The K+/Na+-exchanger uses energy (ATPase) to exchange 2 potassium ions with 3 sodium ions. We 
used a pump-rate per unit area defined by: 

1.5 

The maximum pump rate, Pmax^ and the threshold concentrations, KNai and Kke^ are given in Table 
1. 

Initial conditions 

Initial conditions were determined in the following way: As a starting point, we used c/eno = CknL and 
'^Mo = vml as our initial conditions, where CknL and vml were the resting concentrations and resting 
membrane potential found in a previous study [15]. We then ran a simulation with no system input or 
output. With the membrane mechanisms included in Eqs.[3T]-[33l the system had a simulated resting state 
{cknS and vms) which was close to, but not identical to CknL and vml- For all subsequent simulations, 
we set the initial conditions to the simulated resting conditions (c/eno = CknS and vmq = vms)- The 
estimated values and the values from the literature are given in Table 1. 
Prior to all simulations, we defined the static charge densities: 

PsI = CmVmO - F{ckIO + CNalO - Cciio) (37) 

psE = CmVmO " F{ckEO + CNaEO " CciEo)- (38) 

aE 
The static charge densities ensure that the total charge density in / and E are consistent with vmo^ 
according to Eq. [71 

Comparison of concentrations and charges 

To allow direct comparison with ion concentrations, we represented the charge density in Eq. [71 as an 
equivalent concentration of unit charges, defined by: 

Cen = CKn + CNan " Ccin + psn/F, (39) 

with Eq. [371or Eq. [38lfor p^s- Likewise, we represented the current densities as equivalent unit-charge 
fiux densities, defined by: 

Jen = JKn + JNan ~ Jcin (^^) 

Jen = JKn + JNan ~ Jcin (41) 

Implementation 

The model was implemented in Matlab, and the code will be made publicly available at ModelDB 



(http://senselab.med.yale.edu/modeldb). Simulations were run using the Matlab-solver pdepe^ which 



uses variable time steps. For the simulations presented below, we used a maximum time step of 0.1 s, and 
used 100 segments in the x-direction. Improving the resolution had no visible impact on the predicted 
results. Initial conditions were as listed in Table 1, and the sealed end boundary conditions (Eq. [6]) were 
applied. 
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